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Abstract 



S3 



A practical and new Runge-Kutta numerical scheme for stochas- 
^ " tic differential equations is explored. Numerical examples demonstrate 
the strong convergence of the method. The first order strong conver- 
gence is then proved using Ito integrals for both Ito and Stratonovich 
interpretations. As a straightforward modification of the determinis- 
tic Improved Euler /Heun method, the method is a good entry level 
scheme for stochastic differential equations, especially in conjunction 
with Higham's introduction [SIAM Review, 43:525-546, 2001]. 
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1 Introduce a modified integration 

Nearly twenty years ago Kloeden and Platen [3] described schemes for nu- 
merically integrating stochastic differential equations (sdes). Intervening 
research led to recent developments of useful Runge-Kutta like methods for 
Ito SDEs by Andreas Rossler (TTJ [TO] and for Stratonovich SDEs by Yoshio 
Komori [6j HI [5]. These numerical integration schemes for SDEs are quite 
complicated, and typically do not easily reduce to accurate deterministic 
schemes. This short article introduces a Runge-Kutta scheme for SDEs that 
does straightforwardly reduce to a well known deterministic scheme — the 
variously called Improved Euler, Heun, or Runge-Kutta 2 scheme. 

As well as being a novel practical scheme for the numerical integration of 
SDEs, because of the strong connection to a well known deterministic inte- 
gration scheme, the scheme proposed here serves as an entry level scheme 
for teaching stochastic dynamics. One could use this scheme together with 
Higham's pQ introduction to the numerical simulation of SDEs. Section [2] 
on the method applied to examples assumes a background knowledge of ba- 
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sic numerical methods for ordinary differential equations and deterministic 
calculus as typically taught in early years at university. Section [3] on the 
underlying theory assumes knowledge of stochastic processes such as contin- 
uous time Markov Chains, and, although not essential, preferably at least a 
formal introduction to stochastic differential equations (such as the book [9] 
or article [T] with material that is successfully taught at second/third year 
university) . 

— * 

Consider the vector stochastic process X(t) G M n that satisfies the general 

It6 SDE 

dX = a(t, X) dt + b(t, X) dW, (1) 

where drift a and volatility b are sufficiently smooth functions of their ar- 
guments. The noise is represented by the differential dW which symboli- 
cally denotes infinitesimal increments of the random walk of a Wiener pro- 
cess W(t,u). The symbolic form of the SDE flTJ follows from the most basic 
approximation to an evolving system with noise that over a time step At k 
the change in the dependent variable is 

AX k « a(t k , X k )At k + b{t k , X k )AW k 

where AW k = W(t k+ i,u) — W(t k ,u) symbolises some 'random' effect. This 
basic approximation is low accuracy and needs improving for practical ap- 
plications, but it does form a basis for theory, and it introduces the noise 
process W(t, u), called a Wiener process. We use u to denote the realisation 
of the noise. Such a Wiener process is defined by W(0,u>) = and that the 
increment W(t,u) — W(s,u) is distributed as a zero-mean, normal variable, 
with variance t — s , and independent of earlier times. Consequently, crudely 
put, dW/dt then is a 'white noise' with a flat power spectrum. The SDE (JTJ 
may then be interpreted as a dynamical system affected by white noise. 

The proposed modified Runge-Kutta scheme for the general SDE (jTJ) is the 
following. Given time step h, and given the value X(t k ) = X k , estimate 
X(t k+ i) by X k+ i for time t k+1 = t k + h via 

K x = ha{t k , X k ) + (AW k - S k Vh)b{t k , X k ), 
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K 2 = M(t k+1 , X k + K x ) + (AW k + S k Vh)b(t k+1 , X k + K x ), 

X k+1 = X k + ±(K 1 + K 2 ), (2) 

• where AW k = vhZ k for normal random Z k ~ iV(0, 1); 

• and where S k — ±1 , each alternative chosen with probability 1/2. 

The above describes only one time step. Repeat this time step (t m —t ) / h times 
in order to integrate an sde (Tj[|) from time t — to to t — t m . 

The appeal of the scheme §Z§ as an entry to stochastic integrators is its 
close connection to deterministic integration schemes. When the stochastic 
component vanishes, 6 = 0, the integration step ([2]) is precisely the Improved 
Euler, Heun, or Runge-Kutta 2 scheme that most engineering, science and 
mathematics students learn in undergraduate coursework. 

This connection has another useful consequence in application: for systems 
with small noise we expect that the integration error of the SDE is only a 
little worse than that of the deterministic system. Although Section |3] proves 
the typical 0(ii) error of the stochastic scheme (121), as demonstrated in the 
examples of the next Section |2j when the noise is small expect the error to 
be practically better than the order of error suggests. 

Section[3]also proves that the scheme (j2J) integrates Stratonovich SDEs to 0(ti) 
provided one sets S k = throughout (instead of choosing ±1). 

An outstanding challenge is to generalise this method (T5]) to multiple noise 
sources. 



2 Examples demonstrate O (ti) error is typical 

This section applies the scheme (T5]) to three example SDEs for which, for 
comparison, we know the analytic solution from Kloeden and Platen [3]. 
Two of the examples exhibit errors 0(h), as is typical, whereas the third 
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exhibits a error 0(h 2 }, which occurs for both deterministic odes and a class 
of SDEs. These errors are 'pathwise' errors which means that for any one 
given realisation u of the noise process W(t, u) we refer to the order of error 
as the time step h — > for a fixed realisation co. 



2.1 Autonomous example 

Consider the 'autonomous' SDE 
dX 



iX + Vl + X 2 dt + y/l + X 2 dW with X(0) = 0, (3) 

for some Wiener process W(t, u). The SDE is not strictly autonomous because 
the noise dW introduces time dependence; we use the term 'autonomous' 
to indicate the drift a and volatility b are independent of time. For the 
SDE (jHJ), Kloeden and Platen [3] list the analytic solution as X(t,oS) = 
sinh [t + W(t,u) 



Such analytic solutions are straightforwardly checked via the basic version of 
Ito's formula, 

{df ld 2 f\ df 
if X = f(t, w) for w = W(t, u>), then dX = Uf + --^ dt + dW, 



dt 2 dw 2 J dw 

, (4) 

which may be understood as the usual deterministic derivative rule dX = 
(df /dt) dt + (df /dw) dW with the extra term |(<9 2 / / 'dw 2 ) dt arising from a 
formal multi- variable Taylor series in the infinitesimals dt and dw, recognising 
formally that dW 2 = dt in effect, and all remaining infinitesimal products 
negligible [U EJ e.g.]. 

The proposed numerical scheme (T5]) was applied to integrate the SDE ([3]) from 
t = to end time t = 1 with a time step of h — 1/n for n = 2 16 , 2 15 , . . . , 2 4 
steps. For each of 700 realisations of the noise W(t,u), the Wiener incre- 
ments, AW ~ N(0, 2~ 16 ), were generated on the finest time step, and subse- 
quently aggregated to the corresponding increments for each realisation on 
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Figure 1: As the time step is successively halved, n = 16,32,64,128,256 
time steps over < t < 1 , the numerical solutions of the SDE fl3]) via the 
method (J2J) appear to converge. 



the coarser time steps. Figure [T] plots the predicted X(t,u) obtained from 
the numerical scheme ([2]) for just one realisation u using different time steps. 
The predictions do appear to converge to a well defined stochastic process as 
the step size is repeatedly halved. 

For each size of time step, Figure |5] uses the analytic solution to find the 
RMS error of the predicted X(l,u), averaged over 700 realisations u. This 
RMS error estimates the square-root of the expectation E[(X m — X(1,cj)) 2 ]. 
Figure [2] uses a log-log plot to show that the RMS error decreases linearly 
with time step size h (over four orders of magnitude in time step). That is, 
empirically we see the scheme ([2]) has RMS error 0(fi). 
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-5 

1U _ 5 _ 4 _ 3 _ 2 _1 

10 10 10 10 10 

time step h 

Figure 2: Average over 700 realisations at each of 13 different step sizes for 
the SDE (jHJ): at t = 1, the RMS error in the predicted X(l,co) decreases 
linearly in time step h. 



2.2 Non- autonomous example 



Consider the 'non-autonomous' SDE 



dX 



1 + t 2 



X 2 



dt+(l+t) ( 1 



X 2 



(l + *) s 



3/2 



dW, (5) 



with initial condition that X(0) = , for some Wiener process W(t, u). Here 
both the drift a and the volatility b have explicit time dependence. Ito's 
formula (HI) confirms that the analytic solution to this SDE ([3]) is X(t,oj) = 
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Figure 3: Average over 700 realisations at each of 13 different step sizes for 
the SDE (j5J): at t = 1, the RMS error in the predicted X(l,co) decreases 
linearly in time step h. 



(i + t)w(t, w)/v/iTWp. 

To determine the order of error of the scheme (j2J), the same approach was 
adopted here as described in the previous Section 12.11 The slope of the log- 
log plot in Figure [3] shows that again the RMS error of the predicted X(l,u) 
is Oui) for time step h over four orders of magnitude in h. 
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-l 
10 = 




10 10 10 10 10 



time step h 

Figure 4: Averaging over 700 realisations at each of 13 different step sizes for 
the linear sde (jHJ): at t = 1 , the RMS error in the predicted X(1,cj) decreases 
quadratically, like h 2 . 



2.3 Example with second order error 



Consider the following sde linear in X: 
2X 



dX 



1 + t 



+ (1 + *)' 



dt + (l + tydW withX(0) = l. 



(6) 



for some Wiener process W(t, u). Ito's formula (jH) confirms that the analytic 
solution to this SDE © is X(t, u) = (1 + tf [l + 1 + W(t, u)] . 

To determine the order of error of the scheme (j2J), the same approach was 
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adopted here as described in Section 12.11 The difference is that the slope 
of the log-log plot in Figure H] shows that here the RMS error of the pre- 
dicted X(1,uj) is 0(h 2 ). There appears to be some SDEs for which the error 
of the scheme ([5]) is quadratic in the time step h rather than linear. 

Exercise 1. Use Ito's formula (J4j) to confirm the solutions given below satisfy 
the corresponding given sde. Apply the scheme (j2J) to some of the following 
SDEs and compare the predictions, for different time steps sizes, to the given 
analytic solution. Perhaps adapt some of the code given by Higham [Tj 
Listing 6]. 

1. dX = \(X - t) dt + (X - t - 2) dW, X(0) = 3; solution X = 2 + t + 
exp[W(t)]. 

2. dX = XdW, X(0) = 1; solution X = exp[W(t) - t/2}. 

3. dX = -X(l-X 2 ) dt+(l-X 2 ) dW, X(0) = 0; solution X = tanh[W(t)]. 

4. dX = -Xdt + e-'dW, X(0) = 0; solution X = e-*W(t). 

5. dX = -|X(1 - X 2 fdt + (1 - X 2 fl 2 dW, X(0) = 0; solution X = 

W{t)/y/T+WW- 

For which SDEs is the error 0(h 2 \l 



3 Prove 0(hj global error in general 

This section uses stochastic integration to establish the general order of ac- 
curacy of the proposed numerical integration scheme. 

Proofs that numerical schemes do indeed approximate SDE solutions are often 
complex. My plan here is to elaborate three successively more complicated 
cases, with the aim that you develop a feel for the analysis before it gets 
too complex. Lemma [T] first proves that the Runge-Kutta like scheme (J2J) 
approximates the simplest Ito integrals X = j a b(t) dW to first order in the 
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time step. Second, section I3T21 identifies a class of linear SDEs with additive 
noise when the scheme ([2]) is of second order. Third, section 13.31 proves the 
first order global error of scheme ([2]) when applied to general SDEs. Those 
familiar with stochastic Ito integration could proceed directly to the third 
section 13.31 

One outcome of this section is to precisely 'nail down' the requisite properties 
of the choice of signs Sj in the scheme (j2J). 



3.1 Error for Ito integrals 

This subsection establishes the order of error in computing the Ito integral 
X = £ b(t) dW if one were to invoke the scheme ([2J) on the scalar sde 
dX = b(t) dW. Before proceeding, recall that two fundamental properties 
on the expectation and variance of Ito integrals are widely useful [2J, p. 2] P 
pp. 101-3]: 



martingale property, E 
Ito isometry, E 



f(t,u)dW 

b 

f(t,u) dW 



0; (7) 



E[f(t,oo) 2 ]dt. (8) 



These empower us to quantify errors in the integrals that approximate solu- 
tions of SDEs as in the following lemma. 

Lemma 1. The Runge-Kutta like scheme (j2J) has 0(ti) global error when 
applied to dX = b(t) dW for functions b(t) twice differentiable. 

Proof. Without loss of generality, start with the time step from to — to 
ti = t + h = h . Applied to the very simple SDE dX = b(t) dW one step of 
the scheme d2J) computes 



K x = (AW - SVh)b , K 2 = (AW + sVh)h 
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and then estimates the change in X as 



AX = i(6o + 61) AW + Uh - b )sVh, 



(9) 



where the integrand values 60 — b(0) and b\ = b(h). The classic polynomial 
approximation theorem [7J p. 800, e.g.] relates this estimate (J9]) to the exact 
integral. Here write the integrand as the linear interpolant with remainder: 

bit) = !(&! + b ) + i(6 x - bo)(t - h/2) + \t(t - h)b"(r) 

for some < r(t) < h . Then the exact change in X(t) is 



AX = I b(t) dW = i(6i + 6 )Aiy + i(6 a - 6 ) /' 
Jo Jo 



(t - h/2) dW 



+ 1 



t{t - h)b"(r) dW. 



(10) 



The error in one step of the scheme ([2]) is the difference between the changes ([H 
and ( jTUi) . That is, the true integral change AX = AX + eo where the error 



h - b 
h 



-,Sh 3 / 2 



(t - h/2) dW 



+ ± / t(t-h)b"(r)dW. (11) 



How big is this error? First take expectations, invoke the martingale prop- 
erty (EJ) for the two stochastic integrals, and see that E[e ] = provided 
EfS"] = . Thus the signs S must be chosen with mean zero. 

Second compute the variance of the error e to see the size of the fluctuations 
in the error. Since the expectation E[eo] = 0, the variance Varfeo] = E[eg]. 
Look at various contributions in turn. The first term in the error (ITT]) has 
variance E[(Sh 3 ^ 2 ) 2 } = h 3 E[S 2 ] = 0(h 3 ) provided the signs S have bounded 
variance. Choosing the signs S independently of the noise W there are then 
no correlations between the S terms and the other two terms. The second 
term in the error fill I) has variance 

2 1 rh 

(t - h/2) dW) = J (t- h/2) 2 dt by Ito isometry 



E 
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The third term in the error fill I) , by the Ito isometry (JSJ), has variance 



E 



h ^ 2 

t{t - h)b"(r) dW 



i\t-hfb"(rfdt 



<Blj\\t-hfdt=^Blh\ (12) 



when the second derivative is bounded, < -B2 ■ Lastly, the correlation 

between these previous two integrals is small as, by a slightly more general 
version of the Ito isometry (jSJ), 



E 



(t-h/2)dW / t{t - h)b" (t) dW 



< 



(t - h/2)t(t - h)b"{r) dt 

rh 

B 2 \(t-h/2)t(t-h)\dt = 0(h 4 ). 
Jo 



Hence the local, one step, error is dominated by the first two contributions 
and has variance Var[eo] = 0(h?}. 

To estimate the global integral, J a b(t) dW, we take n — 0(l/h) time steps. 
With n steps the global error is the sum of n local errors: the scheme (J2J) 



approximates the correct solution with global error e = X]j=o e i • Firstly, 
E[e] =0 as Efe^] = for all time steps. Secondly, as the errors on each time 
step are independent, the variance 



Varf 



n-l 

^Var[ ej 

j=0 



n 



Var[e ] = 0(nh 3 ) = 0(h 2 ). 



Thus, for the SDE dX = b(t) dW, the scheme ([2]) has global error of size 0(hj . 



□ 



Tony Roberts, October 4, 2012 



3 Prove 0(h) global error in general 



14 



3.2 Error for linear SDEs with additive noise 



This second lemma addresses somewhat more general scalar SDEs. It not 
only serves as a 'stepping stone' to a full theorem, but illustrates two other 
interesting properties. Firstly, we identify a class of SDEs for which the 
scheme (j2J) is second order accurate in the time step as seen in Example 12.31 
Secondly, the proof suggests that the sign S in the scheme (T5]) relates to sub- 
step properties of the noise W that are independent of the increment AW. 

Lemma 2. The Runge-Kutta like scheme (j2J) has global error 0(h) when 
applied to the additive noise, linear sde dX = a(t)X dt + b(t) dW for func- 
tions a and b twice differentiable. Further, in the exact differential case when 
oh = db/dt (a solution to the SDE is then X = b(t)W) the global error 
isO(h 2 ). 

Proof. In this case, straightforward algebra shows the first step in the scheme (j2J) 
predicts the change 

AX = h\(a + ai)X + |/i 2 a aiX + \{b Q + h)AW 

+ \ ai b Q h(AW - SVh) + \SyfhAb , (13) 

where the coefficient values an = a(0), a\ = a(h), bo = 6(0) and b\ = b{h). 
We compare this approximate change over the time step h with the true 
change using iterated integrals. For simplicity we also use subscripts to 
denote dependence upon 'time' variables t, s and r. Start by writing the SDE 
dX = a t X t dt + b t dW as an integral over the first time step: 

AX = a t X t dt+ / b t dWt 
Jo Jo 

[substituting X t = X + AX inside the first integral] 



h r ft rt 

at X + / a s X s ds+ / b s dW s 
o L Jo Jo 

rh ph pt 

Xq / a t dt + I a t a s X s ds dt 
Jo Jo Jo 



h 



dt+ / b t dW t 



o 
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X + I a r X r dr + b r dW r 
o Jo 



ds dt 



ph pt ph 

+ at b s dW s dt+ / b t dW t 
Jo Jo Jo 

[substituting X s = X + AX inside the second integral] 

ph ph pt 

= X a t dt+ a t a s 
Jo Jo Jo 

ph pt ph 

+ / a t / b s dW s dt+ / btdW t 
Jo Jo Jo 

ph ph pt ph pt ps 

= Xq / a t dt + Xo / at a s ds dt + at a s a r X r dr ds dt 
Jo Jo Jo Jo Jo Jo 

ph pt ps ph pt ph 

+ (k a* K dW r dsdt+ a t b s dW s dt + b t dW t . 
Jo Jo Jo Jo Jo Jo 

(14) 

For the last part of the lemma on the case of higher order error, we need to ex- 
pand to this level of detail in six integrals. Of these six integrals, some signif- 
icantly match the components of the numerical step f[T^]) and some just con- 
tribute to the error. Recall that the proof of Lemma Q] identified that errors 
had both mean and variance. To cater for these two characteristics of errors, 
and with perhaps some abuse of notation, I introduce the notation 0(hP ', h q ) 
to denote quantities with mean 0[hP\ and variance 0(h q ^. For example, 
0(h p , 0) classifies deterministic quantities (9(/i p ), whereas O(0,h q ^) charac- 
terises zero mean stochastic quantities of standard deviation scaling like h q ^ 2 . 
The previous proof looked closely at the variances of error terms; here we 
simplify by focussing only upon their order of magnitude. In particular, 
let's show that the six integrals in ( I14p match the numerical step f lT5|) to an 
error 0(h 3 , h 5 ). 

Consider separately the integrals in ( fl4|) . 

• Firstly, X a t dt = X h^(ao + a±) + C(/i 3 ) by the classic trapezoidal 
rule. This matches the first component in the numerical (1131) . 



Secondly, using the linear interpolation a t = a + + 0{t 2 ), where 
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as usual Aa the double integral 

,dsdt = J (a + (a t + ^t 2 ) + £>(t 3 ) di 

= / alt + aj^t 2 + 0(t 3 )dt 
Jo ^ri 



o Jo 



Aa 



^*h 2 + a ^h 2 + O(h 4 ) 



= la oai h 2 + 0{h A ). 

Multiplied by X , this double integral matches the second term in the 
numerical (ITS]) . 

Thirdly, the triple integral 

r-h rt ps 

/ at a s a r X r dr ds dt = O (h 3 j 
Jo Jo Jo 

because, as seen in the previous two items, each ordinary integration 
over a time of 0(h) multiplies the order of the term by a power of h. 

Fourthly, look at the single stochastic integral in ( Tl4|) . the last term. 
From the proof of the previous lemma, equations fTTUT) and ([121) give 

J b t dW t = \{b 1 + b )AW + — J (t-l)dW t + O(0,h 5 ). (15) 

The first term here matches the third term in the numerical fll3p . The 
second term on the right-hand side is an integral remainder that will 
be dealt with after the next two items. 

Fifthly, change the order of integration in the double integral 
f a t f b s dW s dt 



10 Jo 
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h r h 

b s / a t dtdW s 



o 



I 



b s / ai + 0(h-t)dtdW 8 

/ b s a 1 (h-s) + 0((h-s) 2 )dW s 
Jo 

boa^h-^ + O^dWt 

[by the martingale property d7j) and Ito isometry (|8])] 
h 

b oai (h-t) dW t + O(0,h 5 ) 

|6oai + & ai (| - *) <W t + O (0, /i 5 ) 

= i/i&ociiAiy - 6 oai / + 0(0, /i 5 ) 

Jo 

The first term here matches the first part of the fourth term in the 
numerical f)13p . The second term on the right-hand side is an integral 
remainder that will be dealt with after the last item. 

• Lastly, the triple integral 

/ a t a s b r dW r dsdt = O(0,h 5 ) 
Jo Jo Jo 

because, as in the last item, changing the order of integration to do the 
stochastic integral last, the integral transforms to J (9(/i 2 ) dW which 
by the martingale (0) and Ito isometry (jSJ) is 0(0, h 5 ). 

Hence we now identify that the difference between the Runge-Kutta like 
step (I13p and the change (j!4p in the true solution is the error 

eo = - ±ai& /i 3/2 S + \sVhM) + 6 ai / (* - f ) dW t 

Jo 
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Ab f h 

— — / (t - |) dW t + o(h 3 ) + 0(0, h 5 ) 
h Jo 

= \Sh 3 ' 2 - (t-f)dW t |- a i & o + ^| + 0(h s ,h 5 ). (16) 

Two cases arise corresponding to the main and the provisional parts of 
lemma [2J 

• In the general case, the factor in square brackets, [•], in f)16p determines 
the order of error. Choosing the signs S randomly with mean zero 
then Sh 3 ' 2 = O(0, h 3 ). Recall the integral (t - §) dW t = 0(0, K 3 ) 
also. Thus the leading error is then 0(h 3 ,h 3 y This is the local one 
step error. Summing over 0(l/hj time steps gives that the global 
error is 0(h 2 , /i 2 ). That is, the error due to the noise dominates, vari- 
ance 0(h 2 ), and is generally first order in h as the standard deviation 
of the error is of order h. 

But as the noise decreases to zero, b — > 0, the factor in curly braces, {•}, 
goes to zero. In this decrease the order of error (TIB"]) transitions smoothly 
to the deterministic case of local error 0(h 3 ^j and hence global er- 
ror 0(h 2 ). 



• The second case is when the factor in braces in (flBjl is small: this occurs 
for the integrable case ab = db/dt as then the term in braces is 0(ti) 
so that the whole error (fTB"|) becomes 0(h 3 , /i 5 ). Again this is the local 
one step error. Summing over 0(l/h) time steps gives that the global 
error is 0(h 2 , /i 4 ). That is, in this case the error is of second order in 
time step h, both through the deterministic error and the variance of 
the stochastic errors. Figure H] shows another case when the error is 
second order. 

This concludes the proof. □ 

Interestingly, we would decrease the size of the factor in brackets in the 
error ( TIB"]) by choosing the sign S to cancel as much as possible the integral 
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j (t — |) dW t . This sub-step integral is one characteristic of the sub-step 
structure of the noise, and is independent of AW. If we knew this integral, 
then we could choose the sign S to cause some error cancellation; however, 
generally we do not know the sub-step integral. But this connection between 
the signs S and the integral f Q (t — |) dW t does suggest that the sign S 
relates to sub-step characteristics of the noise process W. 

For example, if one used Brownian bridges to successively refine the numerical 
approximations for smaller and smaller time steps, then it may be preferable 
to construct a Brownian bridge compatible with the signs S used on the 
immediately coarser step sizejj 



3.3 Global error for general SDEs 

The previous section 13.21 established the order of error for a special class of 
linear SDEs. The procedure is to repeatedly substitute integral expressions 
for the unknown whereever it appears (analogous to Picard iteration). In sec- 
tion [372] each substitution increased the number of integrals in the expression 
by two. For general SDEs, this subsection employs the same procedure, but 
now the number of integrals doubles in each substitution. The rapid increase 
in the number of integrals is a major complication, so we only consider the 
integrals necessary to establish that the global error is O (h) . 

Further, the following theorem is also proven for vector SDEs in M. n , whereas 
the previous two subsection sections only considered special scalar SDEs. 

^^The Brownian bridge stochastically interpolates a Wiener process to half-steps in time 
if all one knows is the increment AW over a time step h. The Brownian bridge asserts 
that the change over half the time step, h/2, is ^AIT — \\fhZ for some Z ~ N(0, 1); the 
change over the second half of the time step is correspondingly | AW + \\fhZ . Factoring 
out the half, these sub-steps are \{AW =F ZVh) which match the factors (AW T- S\fh) 
used by the scheme (|2|): the discrete signs S = =f1 have mean zero and variance one just 
like the normally distributed Z of the Brownian bridge. 
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Theorem 3. The Runge-Kutta like numerical scheme (]2]) generally has 
global error 0(h) when applied to the SDE §I§ for sufficiently smooth drift 

— * 

and volatility functions a(t, x) and 6(t, x). 

Proof. The proof has two parts: the first is the well known, standard, ex- 
pansion of the solution of the general SDE (TjQ) by iterated stochastic inte- 
grals leading to the Milstein scheme [U EJ e.g.]; the second shows how the 
scheme ([2]) matches the integrals to an order of error. 

First look at the repeated integrals for one time step; without loss of gener- 
ality, start with a time step from to — to t\ = to + h = h as the analysis 
for all other time steps is identical with minor shifts in the times of eval- 
uation and integration. The stochastic Taylor series' analysis starts from 
the integral form of Ito formula (jl]): for a stochastic process X(t) satisfying 
the general Ito SDE (TjQ), for operators L*, any smooth function f(t,x) of the 
process satisfies 

f(t,X t ) = f(0,j£ o )+ f L°J(s,X s )ds+ f Llf(s,X 8 )dW 8 , (17) 

Jo Jo 



where L 



d_ 
dt 







d 2 



1 dx 



1 j t 



For conciseness we use subscripts 0, t, s and r to denote evaluation at these 
times, and similarly f t = f(t,X t ), and use subscripts % and j to denote 
components of a vector, with the Einstein summation convention for repeated 
indices. As you would expect, when stochastic effects are absent, 6 = 0, the 
integral formula ( flTI) reduces, through the first two components of L°, to 
an integral version of the well known deterministic chain rule: f(t,X t ) = 
f{0, X ) + j l [d t f(s, X s ) + aid x J{s, X s )] ds . Now turn to the SDE (HD itself: 
it is a differential version of an integral equation which over the first time 
step gives 



AX = X(h,u)-X(0,u) 



dX 
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a(t, X t ) dt 



b(t,X t ) dW t 



[apply the Ito formula (IT7|) to both a(t, X t ) and b(t,X t )] 



ph r ft ft 

j a + L° s a s ds + L\a s dW s 
Jo L Jo Jo 



dt 



+ 



h r 



6 + / L%ds+ / L&dW 



[apply the Ito formula i fPTj) to L^fe, 



ft 



+ 



+ 



c?o dt 
ft 




L°a s ds d£ + 



ft, 



o ./o 
b dW t + 




L]a, dW q dt 



o Jo 



h ft 



o 




L°6, ds dW t 



o ./o 




JO 



L$o 



L° r Llb r dr + 



L\L% dW r 



dW, dW+ 



[now rearrange these eight integrals in order of magnitude] 



a 
+ 



fh fh ft 

dt + b / dW t + L^bo / / dW s dWt 
Jo Jo Jo 



h ft 




Lla, dW. dt 



'0 Jo 

+ 




L% dsdW t 




+ 



JO JO 
h ft 






ft 

JO 

L\L% dW r dW s dW t 

h ft 



LM S ds dt + 



o Jo 




L° r Llb r dr dW s dW t 



JO JO 



(18) 



Simplify the first line in this last expression (|18p for AX using the 
well known integrals f Q dt = h, f dW t = AW and f Q f* dW s dW t = 

J^W t dW t = \(AW 2 - h) [1, (3.6), e.g.]. The last of these three in- 
tegrals follow from applying Ito's formula applied to F(t, W t ) = \W 2 
to deduce dF = \ dt + W t dW t , and integrating a rearrangement gives 
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JW t dW t = J dF - j\dt = \W 2 -\t. Also simplify the first line by 
defining the matrix b' Q = [dbi/ dxj\ t _ so that L\bo = b' bo ■ 

• The three integrals above in square brackets in expression (|18p all have 
expectation zero and variance (9(/i 3 ). Recall that with two arguments 
0(h p ,h q ) denotes quantities with mean 0(h p ) and variance 0(h q ^). 
Thus these three integrals in square brackets are O(0, h 3 ). 

• The two integrals above in curly braces in expression (fT8l) are all 0(h 2 } 
in magnitude and hence are 0(h 2 , /i 4 ). 

Combining all these leads to the well established Milstein scheme for the 
change in X over one time step from to to t\ as 

AX = d h + b AW + b' b ±{AW 2 -h) + 0(h 2 , h 3 ) . (19) 



Second, we proceed to show the scheme (jzj) matches this Milstien scheme (1191) . 
Note K x = ha + (AW-SVh)b = 0(h, h) so the product K X K X = 0(h, h 2 ) 
and so on. Hence, by Taylor series in the arguments of the smooth drift a 
and volatility 6, 



K 2 = h a + d!^K x + 0(h, h 2 ) 

+ (AW + SVh) \b + hb + b'oKx + \%K y K x + 0(h 2 , h 3 ) 



where b"K\K\ denotes the tensorial double sum d 2 b/dxidxjKuKij, and 

where the overdot denotes the partial derivative with respect to time, b = 
db/dt. Combining K\ and K 2 , the corresponding first step in the scheme d2J) 
predicts the change 

AX = d h + boAW + \b' b (AW 2 - S 2 h) 

+ \{AW - SVh) \hd' bo + \(AW 2 - S 2 h)b'$ bo 

+ \h(AW + SVh)(b + b' d ) + 0(h 2 , h A ) . (20) 
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Provided S 2 = 1 + 0[h,h) the first lines match to 0(h 2 ,h 3 y. normally 
S 2 = 1 as specified in fl2j). Other terms detailed in ( |2"Uj) are O\0, h 3 ) provided 
E(S') = 0(0, h): normally set to be zero as specified in ([2]). Hence one step 
of the scheme ([2]) matches the solution to 0(h 2 , /i 3 ) . The local error over one 
step of 0(h 2 , h 3 ) leads to, over 0{l/li) steps, a global error of 0(h, h 2 ). □ 

This proof confirms the order of error seen in the earlier examples. Further, 
because we can readily transform between Ito and Stratonovich SDEs, we now 
prove that a minor variation of the numerical scheme applies to Stratonovich 
SDEs. 

Corollary 4 (Stratonovich SDEs). The Runge-Kutta like scheme (|2]) ; but 
setting S = , has errors 0(h)j when the SDE (CQ) is to be interpreted in the 
Stratonovich sense. 

Proof. Interpreting the SDE ([T]) in the Stratonovich sense implies solutions 
are the same as the solutions of the Ito SDE 

dX = {a + \b'b)dt + bdW. 

Apply the scheme ([2]) (with S = ±1 as appropriate to an Ito SDE), or the 
analysis of the previous proof, to this Ito SDE. Then, for example, the one 
step change (12U|) becomes 

AX = (a + %b Q )h + b AW + lb' b (AW 2 -h) + 0(h 2 , h 3 ). 

The component of the deterministic drift term that involves bob' cancel leav- 
ing, in terms of the coefficient functions of the Stratonovich SDE, 

AX = a h + b AW + \b$oAW 2 + 0(h 2 , h 3 ). (21) 

Now apply the scheme ([2]) with S — to the Stratonovich SDE: Taylor series 
expansions obtain the one step numerical prediction as fl20|) upon setting 
S = . This one step numerical prediction is the same as (|2"T|) to the same 
order of errors. Thus the scheme §Z§ with 5* = solves the Stratonovich 
interpretation of the SDE ([T]). □ 
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Exercise 2 (iterated integrals). Consider the scalar sde dX = X dW . This 
SDE is shorthand for the Ito integral X t = X + JqX s dW s . Over a small time 

interval At = h this integral gives Xh = Xq + J^X t dWt ■ Use this as the 
start of an iteration to provide successively more accurate approximations 
to Xh- successive approximations are successive truncations of 

X h « X + X AW + X [\{AWf - \h] + X [|(AW) 3 - \hAW] . 

Determine the integral remainders for each of the approximations. 

Exercise 3 (quadratic convergence). Adapt the proof of Lemma [2] to prove 
that in the specific case when the drift a = a(t) + fi(t)X and the volatility, 

independent of x, satisfies b — (3b, then the scheme has local error 0(h 3 , h 5 ) 
and hence global error O^h 2 , /i 4 ), as seen in Figure IU 

4 Conclusion 

A good basic numerical scheme for integrating Ito SDEs is the Runge-Kutta 
like scheme 02]) (set Sk = to integrate Stratonovich SDEs). A teacher could 
introduce it in the context of the introduction to numerical SDEs outlined by 
Higham pp. 

One of the appealing features of the scheme (j2D is that it reduces, for small 
noise, to a well known scheme for deterministic odes. Consequently, we 
expect the global error (9(||a||/i 2 + \\b\\ti) for some norms of the drift and 
volatility. Such more general expressions of the error should be useful in 
multiscale simulations where the strength of the noise depends upon the 
macroscale time step, such as in the modelling of a stochastic Hopf bifurca- 
tion [8, §5.4.2]. 

One required extension of the scheme fl2]) is to generalise it, if possible, to 
the case of multiple independent noises. I am not aware of an attractive 
generalisation to this practically important case. 
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